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Discussion of "Statistical Modeling of 
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We congratulate the authors for producing such 
a helpful and comprehensive overview paper of a ra- 
pidly developing and important area. The starting 
point for inference in spatial extreme value prob- 
lems is to identify which features of the process re- 
quire modeling. In certain applications, for exam- 
ple, in the generation of return period maps, only 
the marginal behavior is of concern. For such ap- 
plications, the covariate hierarchical/latent variable 
models reviewed in Section 4 are ideal. However, 
if there is any dependence in the process between 
sites, then care needs to be taken when assessing 
the uncertainty in the estimated marginal distribu- 
tion; the composite likelihood procedures detailed 
in Section 6.2 can also be exploited in this context 
when one wishes to avoid assumptions on the form 
of the spatial dependence. As the authors point out, 
however, if interest lies in modeling the joint oc- 
currence of extremes over a region, then the depen- 
dence structure of the spatial process needs to be 
explicitly modeled. The most widely used approach 
in such cases is to model the process as a max- 
stable process. Here we will explore the suitabil- 
ity of this framework for modeling spatial extremes, 
since these processes are quite restrictive in their 
assumptions. Specifically, all finite dimensional dis- 
tributions of a max-stable process are multivariate 
extreme distributions. Even simply in the bivariate 
case, this corresponds to the variables being exactly 
independent or asymptotically dependent. Conse- 
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quently, the broad class of asymptotically indepen- 
dent variables is precluded under the modeling as- 
sumptions of max-stable processes. 

First consider diagnostic testing for the process 
being max-stable. From our experience we feel that 
in many applications, max-stable processes are as- 
sumed to be appropriate without testing their suit- 
ability for the data. A partial justification for this 
is that often it is pointwise maxima that are being 
modeled, and so just as one appeals to the marginal 
limit theory to justify fitting the GEV distribution 
site-wise, it seems natural to appeal to the limit 
theory for the dependence structure also. However, 
we would typically not fit the GEV to the margins 
blindly, but look to assess its suitability through di- 
agnostics such as Q-Q plots. To our knowledge, there 
currently are no diagnostics for testing if the process 
is max-stable, and so this discussion contribution 
aims to provide such a test. 

Suppose that {Y{x) : x £ A} is a max-stable pro- 
cess on the region A with standard Gumbel marginal 
distributions for all x. This is simply attained through 
a pointwise log transformation of the more commonly- 
assumed standard Frechet margins. As the process 
will typically only be observed at a finite set of 



sites xi, . 



then all that can be tested for in 



practice is that the joint distribution of {Y(xj) : 
j € A}, where A = {1, . . . , m}, follows a multivariate 
extreme value distribution. Then for any DCAwe 
have that the joint distribution function for {Y(xj) : 
j G D} is 

Gd(y) = exp[-Vb{exp(y)}], 

where Vd is the associated exponent measure; see 
Section 2.3. A key property of Gd, due to max-sta- 
bility, is that the distribution of Yd = maxjgD y(x,) 
is 

Hd(v) = G D (yl) =exp[-exp{-(y-/io)}]. 

This is a Gumbel distribution with location param- 
eter fi£) = logVb(l) where < fin < log(|D|), due 
to bounds on the exponent measure. It follows that 
Zd = Yd — fj,£> is standard Gumbel for all DC A. 
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The idea of the diagnostic is to pool values of Zry 
over replicates of the max-stable process and over 
all D C A, of a particular cardinality k = \D\, and 
test using a P-P plot whether the variables Zjj fol- 
low a standard Gumbel distribution. This enables 
an assessment of the kth. dimensional properties of 
the process. Here we look at k = 2, 3, 4, m. 

There are a few practical issues to address. As 
the value of /i£> is unknown for all D with \D\ > 2, 
these parameters require estimation. We used max- 
imum likelihood to estimate [id, based upon repli- 
cate data for Yd, and subject to the parameter con- 
straints above; this was found to give better esti- 
mates than moment methods. There are (^|) sets 

with size |D| in the power set of A, thus, for m 
large, it is computationally intensive to determine 
an estimate for all possible /id- As a result, we limit 
ourselves to a randomly selected sample of 500 pos- 
sibilities for D with |L>| = k. Although we did not 
explore this choice in more detail, it would seem 
reasonable not to examine subsets D in which the 
included sites are likely to be independent, since the 
ability to detect a departure from max-stability will 
be limited in these cases. In order to determine the 
variability of our estimates, we use a nonparamet- 
ric bootstrap to produce 95% confidence intervals. 
Bootstrap methods are required as there is clearly 
dependence in the pooled data over different D. By 
treating each replicate of the process as a block, and 
constructing Zd through estimation of ftp for each 
bootstrap sample, this complicated dependence is 
naturally incorporated into the confidence intervals. 

Three different simulated data sets were used to 
illustrate the methods, with 1000 replicates of the 
variables generated on a regular grid of sites over 
a 10 x 10 unit square. The different dependence struc- 
tures used were as follows: a Smith model max- 
stable process (with identity covariance matrix il); 
a multivariate extreme value distribution with lo- 
gistic dependence structure (dependence parameter 
a = 0.7); and a Gaussian process with exponential 
correlation function (A = 0.7" 1 ), pointwise trans- 
formed to have standard Gumbel marginal distri- 
butions. The first two examples are max-stable in 
their dependence structure, that is, follow multivari- 
ate extreme value distributions at the sites on the 
grid. The third example, having a Gaussian copula, 
is not max-stable. 

Figures 1-3 show the diagnostic P-P plots for the 
three dependence structures respectively. Each fig- 
ure illustrates the diagnostic for k = 2, 3, 4, 100, with 
the P-P plot presented on the y-axis as a difference 





Fig. 1. Rescaled P-P plot for Zd derived from the Smith 
max-stable process model. Top-bottom: \D\ =2,3,4,100. 



Fig. 2. As Figure 1, but for the multivariate logistic distri- 
bution. 

between the empirical and model probabilities. For 
the first two processes the horizontal line, represent- 
ing no departure from the model probability, falls 
inside the pointwise 95% confidence intervals, thus 
indicating that the data provide no evidence against 
a max-stable model. This is not the case for the 
Gaussian process copula, which produces evidence 



DISCUSSION 
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Fig. 3. As Figure 1, but for the Gaussian copula. 

to reject max-stability for each value of k. Interest- 
ingly, the power of the test increases with k (note 
the change of scale on the bottom subplot). If only 
the k = 2 case had been used with a smaller sample 
size or weaker dependence we would wrongly have 
failed to reject the null for the Gaussian copula. 
This indicates that higher order dependence plays 
an important role in identifying the nature of the 
extremal dependence. Currently, however, the stan- 
dard methods to fit max-stable processes use only 
the bivariate distributions, through composite like- 
lihood methods, as detailed in Section 6.2. 

More generally, this highlights that pairwise mod- 
eling techniques may not always be as effective as 
one would hope. Even if the process has bivariate ex- 
treme value pairwise distributions, this does not en- 
sure that higher order distributions are max-stable, 
or even asymptotically dependent. The restriction 
to pairwise likelihood is to a large part due to in- 
tractability of higher order distributions of max- 
stable processes. As noted by the authors, higher 
order distributions for the Smith model can be ex- 
pressed (Genton, Ma and Sang (2011)), though the 
model's utility is limited due to its unrealistic pro- 
cess realizations. In terms of the copula models of 
Section 5, Nikoloulopoulos, Joe and Li (2009) pro- 
vide the (i-dimensional copulas for the extremal t 
and, through the appropriate limits, the Hiisler- 
Reiss analogue. The authors comment on the in- 
triguing links between the copula and process mod- 
els. It would seem to us that the extremal t cop- 



ula has the immediate interpretation as the limit- 
ing finite-dimensional dependence structure of nor- 
malized spatial t processes. Kabluchko, Schlather 
and de Haan (2009) similarly demonstrate that the 
Brown-Resnick, analogous to the Hiisler-Reiss, limit 
arises from normalized maxima of Gaussian pro- 
cesses which become increasingly dependent through 
contraction of the spatial domain at an appropriate 
rate. Consequently, these two copula models could 
equally be viewed as max-stable process models. We 
thus agree with the authors' suggestion that the con- 
nection is indeed a matter of extending the copula 
to the full spatial domain. Evidently, the utility of 
full spatial process models lies in being able to as- 
sess features at unobserved sites; this is aided in 
the current context by simulation of the full pro- 
cess. For the extremal t model, the definition-based 
method described in Section 7.1 and illustrated in 
Figure 7 seems adequate for this. The alternative 
would be to derive the appropriate spectral process 
W(x) and exploit representation (20), though it is 
not clear that this would be simpler. For the Hiisler- 
Reiss extension, it appears that correspondence with 
Brown-Resnick processes would permit simulation, 
through the methods already employed in the paper. 

At a more conceptual level, the major restriction 
of max-stable processes is that the form of the ex- 
tremal dependence at observed levels must be as- 
sumed to hold for all more extreme events. How- 
ever, it may often seem more plausible that the de- 
pendence could weaken at extreme levels, with the 
largest extremes becoming more isolated as energy 
is concentrated into a more localized area. While 
max-stable process models may well be sufficiently 
flexible to describe observed extremal dependence, 
for the purposes of extrapolation we would really like 
to be sure that the assumed stability holds. This can 
also be examined through the empirical decay of tail 
probabilities. Suppose X(x) represents the original 
data process (before taking pointwise maxima), but 
transformed to have common Gumbel margins. De- 
fine -Fd(xI) = Pr{X(xj) > x:i € D}; if the depen- 
dence structure is asymptotically dependent, then 
Frj(xl) = 0{exp(— x)} as x — > oo. If this feature is 
observed in the data, this suggests that the limit- 
ing max-stable process is different from indepen- 
dence, and extrapolations based upon fitting these 
processes may be reliable. If asymptotic indepen- 
dence is present in the data, Fd(x1) = 0{exp(— xj 
rjo)}, for some t\b 6 (0, 1). If this feature is observed 
in the data, then we would be rightfully skeptical that 
fitting a max-stable process would provide reliable 
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extrapolations. Some recent work by Wadsworth and 
Tawn (2012) explores the extremal properties of a 
class of spatial processes which satisfy this property. 
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